library(plotly)
library(Hmisc)
library(WVPlots)
library(plyr)
library(reshape)
source("tools/DataSplitTools.R")
source("tools/GeneralizedNadarayaWatson.R")
source("tools/CommonTools.R")
source("tools/CrossValTools.R")

Аналіз даних ЗНО 2018 року та статистики рідних мов

Зчитаємо дані:

df <- read.csv2(file = "data/ZNOandVoating/input_2018.csv", header = TRUE, fileEncoding = "CP1251", sep = ",", dec=".", row.names = 1)
names(df) <- c("region", "ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")
df <- df[, -(4:8)]

describe(df)
## df 
## 
##  3  Variables      315442  Observations
## ---------------------------------------------------------------------------
## region 
##        n  missing distinct 
##   315442        0       23 
## 
## lowest : Вінницька область        Волинська область        Дніпропетровська область Житомирська область      Закарпатська область    
## highest: Херсонська область       Хмельницька область      Черкаська область        Чернівецька область      Чернігівська область    
## ---------------------------------------------------------------------------
## ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   315442        0       84    0.997    117.5    62.76        0        0 
##      .25      .50      .75      .90      .95 
##      105      131      159      180      187 
## 
## lowest :  -1.0   0.0 100.0 102.0 104.0, highest: 197.5 198.0 199.0 199.5 200.0
## ---------------------------------------------------------------------------
## math 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   315442        0       55    0.675    35.38     56.2       -1       -1 
##      .25      .50      .75      .90      .95 
##       -1       -1      100      150      170 
## 
## lowest :  -1   0 100 104 107, highest: 196 197 198 199 200
## ---------------------------------------------------------------------------

Змерджимо із інформацією про рідну мову по всім областям. Отримаємо:

mt <- read.csv2("data/mother_tongue/grouped_mt.csv", header = TRUE, sep = ",", dec=".", row.names = 1)
mt$region <- rownames(mt)
mt$region[mt$region == "м. Київ"] <- "м.Київ"

colnames(mt)[1] <- "ukr_lang"
df[df$math < 0, c("math")] <- 0

df <- merge(df, mt, by = "region", all.x = TRUE)

rownames(df) <- 1:nrow(df)

ЗНО з математики та української мови

Виведемо загальні характеристики даних ЗНО з математики та української мови:

describe(df[, c("math", "ukr")])
## df[, c("math", "ukr")] 
## 
##  2  Variables      315442  Observations
## ---------------------------------------------------------------------------
## math 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   315442        0       54    0.588    36.07    55.77        0        0 
##      .25      .50      .75      .90      .95 
##        0        0      100      150      170 
## 
## lowest :   0 100 104 107 111, highest: 196 197 198 199 200
## ---------------------------------------------------------------------------
## ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   315442        0       84    0.997    117.5    62.76        0        0 
##      .25      .50      .75      .90      .95 
##      105      131      159      180      187 
## 
## lowest :  -1.0   0.0 100.0 102.0 104.0, highest: 197.5 198.0 199.0 199.5 200.0
## ---------------------------------------------------------------------------

Із таблиці вище видно, що в записах нема пропущених значень. Судячи із значень квантилів, незначна частка результатів ЗНО з математики є нулями (менше 10%). Варто відмітити, що трішки більша частка значень 0 у ЗНО з української мови (бфльше 10%, але менше 25%).

Подивимося на розподіли:

p1 <- plot_ly(x = df[, "math"], type = "histogram", name = "math")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram", name = "ukr") 

plotly::subplot(p1, p2) %>% layout(title = "EIT: Ukranian language and mathematics")
remove(p1); remove(p2)

Розподіли балів результатів ЗНО з української та математики помітно різняться.

Подивимося на співвідношення людей, котрі отримали 0 з математики та українській одночасно.

print("Не подолали поріг або не здавали ЗНО відносно всіх абітурієнтів: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df), 2)) %&% "%")
## [1] "Не подолали поріг або не здавали ЗНО відносно всіх абітурієнтів: 13.27%"
print("Не подолали поріг або не здавали ЗНО відносно тих, хто не склав математику: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df[(df$math == 0), ]), 2)) %&% "%")
## [1] "Не подолали поріг або не здавали ЗНО відносно тих, хто не склав математику: 17.83%"
print("Не подолали поріг або не здавали ЗНО відносно тих, хто не склав українську: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df[(df$ukr == 0), ]), 2)) %&% "%")
## [1] "Не подолали поріг або не здавали ЗНО відносно тих, хто не склав українську: 96.21%"

Із виведеної аналітики напрошуються наступні висновки, 1) Невелика частка абітурієнтів, що мають 0 одночасно з двох предметів. 2) Абітурієнти, котрі не склали іспит із математики, скоріш за все не склали й українську.

Видалимо стрічки, що відповідають спостереженням із 0 хочаб в одному ЗНО та подивимося на розподіл.

df <- df[(df$math != 0) & (df$ukr != 0) ,]
p1 <- plot_ly(x = df[, "math"], type = "histogram", name = "math")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram", name = "ukr")

plotly::subplot(p1, p2) %>% layout(title = "Distribution of the passed EIT",
  annotations = list(list(x = 0.2, y = 1.05, text = "math", showarrow = F, xref = 'paper', yref = 'paper'),
                     list(x = 0.8, y = 1.05, text = "ukr", showarrow = F, xref = 'paper', yref = 'paper')))
remove(p1); remove(p2);

Відмітимо, що наразі частка людей, котрі погано здали українську мову зменшилася близько вдвічі (було в околі 900, стало в околі 400).

Гістограма розсіювання, в свою чергу, виглядає наступним чином:

plot_ly(data = df, x = ~math, y = ~ukr, mode = "markers", type = "scatter",  marker=list(opacity = 0.7)) %>% layout(title = "Passed EIT results")

Варто відмітити, що для абітурієнтів, що склали математику краще ніж 190 балів, помітна тенденція, як правило, кращого складання ЗНО з української.

Рідна мова

Виведемо загальні характеристики кожної політичної сили:

describe(df[, -c(1, 2, 3)])
## df[, -c(1, 2, 3)] 
## 
##  5  Variables      78987  Observations
## ---------------------------------------------------------------------------
## ukr_lang 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    78987        0       23    0.996    77.69    19.54    46.28    50.20 
##      .25      .50      .75      .90      .95 
##    67.00    81.00    94.81    97.00    97.26 
## 
## lowest : 46.28362 50.20360 53.80113 67.00188 69.19777
## highest: 95.32371 97.00455 97.26233 97.80575 98.33764
## ---------------------------------------------------------------------------
## rus 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    78987        0       23    0.996     19.8    17.83    2.506    2.729 
##      .25      .50      .75      .90      .95 
##    4.735   15.504   31.914   44.294   44.294 
## 
## lowest :  1.191217  1.776793  2.506493  2.728681  2.902247
## highest: 29.257209 31.914224 41.947480 44.293571 48.194477
## ---------------------------------------------------------------------------
## mold 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    78987        0       23    0.996   0.4461   0.7547  0.01623  0.01623 
##      .25      .50      .75      .90      .95 
##  0.02548  0.03882  0.08728  0.57930  3.77657 
##                                                                       
## Value       0.02  0.03  0.04  0.05  0.06  0.09  0.10  0.15  0.44  0.58
## Frequency  15288 21229  7029 13136  2448  3130  3280  2137  1770  2359
## Proportion 0.194 0.269 0.089 0.166 0.031 0.040 0.042 0.027 0.022 0.030
##                       
## Value       3.78  6.79
## Frequency   6141  1040
## Proportion 0.078 0.013
## ---------------------------------------------------------------------------
## bilorus 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    78987        0       23    0.996   0.1145  0.04561  0.04758  0.05505 
##      .25      .50      .75      .90      .95 
##  0.09161  0.11686  0.13279  0.17519  0.17541 
## 
## lowest : 0.03495828 0.04544391 0.04758436 0.04929121 0.05504710
## highest: 0.15998877 0.16534317 0.17519257 0.17540883 0.19230380
## ---------------------------------------------------------------------------
## other 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##    78987        0       23    0.996    1.947    2.389   0.1332   0.3058 
##      .25      .50      .75      .90      .95 
##   0.3526   0.8635   1.7601   7.8832   7.8832 
##                                                                       
## Value       0.10  0.15  0.25  0.30  0.35  0.40  0.50  0.55  0.70  0.80
## Frequency   2357  2820  2609  3280  8743  1671  4303  2504  2448  5907
## Proportion 0.030 0.036 0.033 0.042 0.111 0.021 0.054 0.032 0.031 0.075
##                                                                 
## Value       0.85  1.05  1.45  1.60  1.75  2.45  7.90 12.30 16.00
## Frequency  10382  2533  3626  2137  6692  8280  6141  1040  1514
## Proportion 0.131 0.032 0.046 0.027 0.085 0.105 0.078 0.013 0.019
## ---------------------------------------------------------------------------
colors_ <- c("aliceblue", "antiquewhite", "aqua", "aquamarine", "azure",
             "beige", "bisque", "black", "blanchedalmond", "blue",
             "blueviolet", "brown", "burlywood", "cadetblue",
             "chartreuse", "chocolate", "coral", "cornflowerblue",
             "cornsilk", "crimson", "cyan", "darkblue", "darkcyan",
             "darkgoldenrod", "darkgray", "darkgrey", "darkgreen")


pivot_mt <- mt[, -ncol(mt)]

n_col <- ncol(pivot_mt)
n_row <- nrow(pivot_mt)

all_plots <- list()
for (i in 1:n_row) {
  all_plots[[i]] <- plot_ly(pivot_mt, 
                            x = colnames(pivot_mt)[-1], y = as.numeric(pivot_mt[i, 2:n_col]), 
                            type = "bar", name = pivot_mt[i, 1], marker = list(color = colors_[1:(n_col-1)]))
}

plotly::subplot(all_plots, nrows = 4, shareX = TRUE) %>% 
  layout(title="Distributions according to regions")
all_plots <- list()
all_plots[[1]] <- plot_ly(x = rownames(pivot_mt), y = pivot_mt[, 1], 
                            type="bar", name = colnames(pivot_mt)[1])
for (i in 2:n_col) {
  all_plots[[i]] <- plot_ly(x = rownames(pivot_mt), y = pivot_mt[, i], 
                            type="bar", name = colnames(pivot_mt)[i])
}

plotly::subplot(all_plots, nrows = 4, shareX = TRUE) %>% 
  layout(title="Distributions of every language")

Розбиття

Розіб’ємо вибірку на тренувальну та тестову частини у відсотковому співвідношенні 80/20%. Далі, тренувальну частину розіб’ємо на 5 рівних частин для проведення 5-ти фолдової кроссвалідації, аналогічно до того, як було зроблено для результатів ЗНО 2016 року.

set.seed(42)

df[, -c(1, 2, 3)] <- df[, -c(1, 2, 3)]/100
df <- df[, -1]
df <- df[(df$math != 0) & (df$ukr != 0) ,]

train_test_list <- train_test_split(df, ratio = 0.80)

train <- train_test_list[["train"]]; test <- train_test_list[["test"]]
write.csv(train, "data/computation_results/train2018.csv"); write.csv(test, "data/computation_results/test2018.csv")

remove(df); remove(train_test_list); remove(test)

Узагальнені оцінки Надарая-Ватсона

Визначення параметрів згладжування

Скористаємося 5-ти фолдовою кросс-валідацією для знаходження оптимальних параметрів згладжування для кожної моделі аналогычно до того, як це було зроблено для результатыв ЗНО 2016 року.

# h_range <- c(seq(0.1, 1, 0.5), seq(1, 4, 0.5),5,  8, 10)
# 
# for (i in 1:5) {
#   print("***** " %&% as.character(i) %&% " ******")
#   set.seed(i)
#   cv_split <- cross_validation_split(train, k = 5)
#   GNW_cv_results <- GNWcv_across_h(h_range = h_range,
#                                  cv_df_split = cv_split,
#                                  X_colname = "math", Y_colname = "ukr",
#                                  W_colname = colnames(cv_split[[1]])[-(1:2)],
#                                  use_parallel = TRUE)
#   write.csv(GNW_cv_results, "data/computation_results/2018mt_data/5cv_mean_std_" %&% as.character(i) %&% ".csv")
# }
# 
# # h_opt <- optimal_h(GNW_cv_results)
# remove(h_range); remove(i); remove(cv_split); remove(GNW_cv_results)

Результат середніх та диспресії збережемно в змінній \(GNW\_cv\_results\), а оптимальні \(h\), серед запропонованих для перебору, у змінну \(h\_opt\).

data_path <- "data/computation_results/2018mt_data/"
file_names <- dir(data_path) #where you have your files
file_names <- data_path %&% file_names

cv_results <- lapply(file_names, 
                     function(x){
                       df <- as.matrix(read.csv(x, row.names = 1))
                       l <- list(df=df , h=optimal_h(df))
                       return(l)
                       })

cv_results_h <- lapply(cv_results, function(x){x$h})
cv_results_mean_std <- lapply(cv_results, function(x){x$df})

print(cv_results_h)
## [[1]]
##          [,1]     [,2]     [,3]     [,4]    [,5]
## [1,] 424.1233 437.0176 1978.592 88151.94 2712.17
## [2,]   5.0000  10.0000   10.000    10.00   10.00
## 
## [[2]]
##          [,1]     [,2]     [,3]      [,4]     [,5]
## [1,] 422.8279 425.1279 1440.651 3121136.5 972.0292
## [2,]   2.0000   3.0000    8.000       3.5   5.0000
## 
## [[3]]
##          [,1]    [,2]     [,3]     [,4]     [,5]
## [1,] 422.6143 425.733 1335.551 11013632 3247.552
## [2,]   1.5000   1.500   10.000        8    5.000
## 
## [[4]]
##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 422.8392 438.2502 1300.302 532573.4 872.2502
## [2,]   0.1000  10.0000   10.000      2.5   5.0000
## 
## [[5]]
##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 425.2315 437.0181 1908.857 515915.1 4538.923
## [2,]   8.0000  10.0000    0.100      4.0   10.000

Порівняння результатів кросс-валідації та прогнозування на відкладеній вибірці

Порівняємо тепер результати кросс-валідації з результатми прогнозування на тестовій частині даних.

# train <- read.csv("data/computation_results/train2018.csv", row.names = 1, stringsAsFactors = FALSE)
# test <- read.csv("data/computation_results/test2018.csv", row.names = 1, stringsAsFactors = FALSE)
# h_opt <- cv_results_h[[1]][2, ]
# 
# GNW <- GeneralisedNadarayaWatson$new()
# GNW$train(X_train = as.numeric(train[, "math"]),
#           Y_train = as.numeric(train[, "ukr"]),
#           W_train = as.matrix(train[, -c(1, 2)]))
# GNW$run_cluster()
# prediction_with_acoeff <- GNW$predict_in_parallel(X_test = as.numeric(test[, "math"]),
#                                                   W_test = as.matrix(test[, -c(1, 2)]), h = h_opt)
# GNW$stop_cluster()
# 
# prediction <- prediction_with_acoeff[["prediction"]]
# A_coeff <- prediction_with_acoeff[["A_test"]]
# 
# write.csv(prediction, "data/computation_results/prediction2018m.csv")
# write.csv(A_coeff, "data/computation_results/A_coeff2018m.csv")
# 
# remove(GNW); remove(train)

Порахуємо зважену суму МНК для прогнозу зробленого на тестових даних.

train <- read.csv("data/computation_results/train2018.csv", row.names = 1, stringsAsFactors = FALSE)
test <- read.csv("data/computation_results/test2018.csv", row.names = 1, stringsAsFactors = FALSE)

prediction <- as.matrix(read.csv("data/computation_results/prediction2018m.csv", row.names = 1, stringsAsFactors = FALSE))
A_coeff <- as.matrix(read.csv("data/computation_results/A_coeff2018m.csv", row.names = 1, stringsAsFactors = FALSE))

WMSE <- weighted_MSE(Y_true = as.numeric(test[, "ukr"]), Y_predicted = prediction, A_coeff = A_coeff)
names(WMSE) <- rep("comp.") %&% as.character(1:5)
WMSE
##      comp.1      comp.2      comp.3      comp.4      comp.5 
##     445.406     421.061   -1484.089 -208261.099   -4124.100

Відмітимо, що найменше значення зваженого МНК для моделей, що відповідають опозиційним силам.

Порівняємо зважені МНК результатів прогнозування з результатами крос-валідації, використовуючи відповідними параметрами згладжування.

h_opt <- cv_results_h[[1]][2, ]
GNW_cv_results <-cv_results_mean_std[[1]]

cv_mean_results <- sapply(1:length(h_opt),
                          function(i){
                            GNW_cv_results[which("mean_" %&% as.character(h_opt[i]) == rownames(GNW_cv_results)), i]
                            })
cv_test_comparing <- rbind(cv_mean_results, WMSE)
rownames(cv_test_comparing) <- c("cross-validation", "test")
colnames(cv_test_comparing) <- 1:5
print(cv_test_comparing)
##                         1        2         3          4        5
## cross-validation 424.1233 437.0176  1978.592   88151.94  2712.17
## test             445.4060 421.0610 -1484.089 -208261.10 -4124.10